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Abstract. We present a new Finite Volume Evolution Galerkin (FVEG) scheme for the 
solution of fhe shallow wafer equations (SWE) wifh fhe boffom fopography as a source 
ferm. Our new scheme will be based on fhe EVEG mefhods presenfed in (Lukacova, 
Noelle and Kraff, }. Comp. Phys. 221, 2007), buf adds fhe possibilify fo handle dry 
boundaries. The mosf imporfanf aspecf is fo preserve fhe posifivify of fhe wafer heighf. 
We presenf a general approach fo ensure fhis for arbifrary finife volume schemes. The 
main idea is fo limif fhe oufgoing fluxes of a cell whenever fhey would creafe negative 
wafer heighf. Physically, fhis corresponds fo fhe absence of fluxes in fhe presence of 
vacuum. Well-balancing is fhen re-esfablished by splitting gravifafional and gravify 
driven parfs of fhe flux. Moreover, a new enfropy fix is infroduced fhaf improves fhe 
reproduction of sonic rarefaction waves. 
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1 Introduction 

The shallow water equations (SWE) are a mathematical model for the movement of water 
under the action of gravity. Mathematically spoken, they form a set of hyperbolic con¬ 
servation laws, which can be extended by source terms like the influence of the bottom 
topography, friction or wind forces. In this case, we will speak of a balance law. For 
simplicity, this work will consider the variation of the bottom as the only source term. 

Many important properties of the model rely on the fact that the water height is 
strictly positive. Despite this, typical relevant problems include the occurrence of dry 
areas, like dam break problems or the run-up of waves at a coast, with tsunamis as the 
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most impressive example. So for simulations of these problems, we have to develop nu¬ 
merical schemes that can handle the (possibly moving) shoreline in a stable and efficient 
way Another crucial point in solving balance laws is the treatment of the source terms. 
For precise solutions, it is necessary to evaluate the source term in such a way that certain 
steady states are kept numerically, i.e. the numerical flux and the numerical source term 
cancel each other exactly for equilibrium solutions. 

In the last years, many groups contributed to the solution of the difficulties described 
above. In [?], Audusse et.al. proposed a reconstruction procedure where the free surface 
and water height are reconstructed and the bottom slopes are computed from these. This 
guarantees the positivity of the water height and gives a well-balanced scheme at the 
same time. Begnudelli and Sanders developed a scheme for triangular meshes including 
scalar transports in [?]. They proposed a strategy how to exactly represent the free surface 
in partially wetted cells, leading to improved results at the wetting/drying front. In [?], 
Brufau et.al. analyse how to deal with flow on an adverse slope. They locally modify the 
bottom topography in certain situations to avoid unphysical run-ups or wave creation at 
the dry boundary. Gallardo et.al. discussed various solutions of the Riemann problem at 
the front and used them in a modified Roe scheme. They then used the local hyperbolic 
harmonic method from Marquina (cf. [?]) in the reconstruction step to achieve higher 
order, see [?]. Kurganov and Petrova proposed a central-upwind scheme that is well- 
balanced and positivity preserving in [?]. It is based on a continuous, piecewise linear 
approximation of the bottom topography and performs the computation in terms of the 
free surface instead of the relative water height to simplify the well-balancing. The last 
feature is also a building block in the work of Liang and Marche [?]. They also provide 
a method to extend this well-balancing feature to situations including wetting/drying 
fronts. Liang and Borthwick [?] used adaptive quad-tree grids to improve the efficiency 
of their schemes. Wetting and drying effects are handled as well as friction terms. In the 
context of residual distribution methods, Ricchiuto and Bollermann developed a positiv¬ 
ity preserving and well-balanced scheme for unstructured triangulations [?]. 

The finite volume evolution Galerkin (FVEG) methods developed by Lukacova, Mor¬ 
ton and Warnecke, cf. [?,?,?], have been successfully applied to the SWF in [?]. They are 
based on the evaluation of so called evolution operators which predict values for the finite 
volume update. Thanks to these operators, the schemes take into account all directions 
of wave propagation, enabling them to precisely catch multidimensional effects even on 
Gartesian grids. These schemes show a very good accuracy even on relatively coarse 
meshes compared to other state of the art schemes and they are also competitive in terms 
of efficiency (cf. [?]). 

However, the existing FVEG schemes are not able to deal with dry boundaries. Thus 
in this work we will present a method to preserve the positivity of the water height with 
an arbitrary finite volume method. To achieve this, we reduce the outflow on draining 
cells such that the water height does not become negative. We will then provide the 
means to preserve the well-balancing property under the presence of dry areas, and apply 
both techniques to a new FVEG method. In addition, we present a new entropy fix for 
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the FVEG schemes that improves the reproduction of sonic rarefaction waves. 

We start our paper with a short presentation of the SWE in Section Section de¬ 
scribes the EVEG method we will start from. The arising difficulties by introducing dry 
areas and means to overcome them are described in Section]^ which is the main part of 
the paper. Einally, in Section we will show selected numerical test cases that demon¬ 
strate the performance of our schemes. 


2 The Shallow Water Equations 

2.1 Balance Law Form 

We consider the shallow water system in balance form 

— -h V-:F(u) = -S (u,x) 


The conserved variables and the flux are given by 

u=lhvX :f(u) = (:Fi(u):p2(u))- 
\hV2 


( hvi 


hV2 


hvl+gY hviV2 
\ hviV2 hvl+g’-Y^ 


( 2 . 1 ) 


( 2 . 2 ) 


where h denotes the relative water height, v= {vi,V 2 y the flow speed and g the (constant) 
gravity acceleration. The source term S (u,x) is given by 


S{u,x)=gh 


db(x) 

db{x) 


(2.3) 


with b{x) the local bottom height. We also introduce the free surface level, or total water 
height, 

H{x) = h{x) + b{x) (2.4) 


and the so-called speed of sound 

c = Vgh- (2.5) 

This is the velocity of the gravity waves and should not be confused with the physical 
sound speed in air. 


2.2 Quasi-linear Form 

For the derivation of the evolution operators in Section [T2j it is helpful to rewrite ( |2.1| | in 
primitive variables. The system then takes the form 

Wf-hAi(w)w;,j-hA2(w)w;c2 = t (2-6) 
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with 



and the source term 


/vi h 0\ 

1 = ^ 0 L ^2 

VO 0 ViJ 



/v2 0 h\ 

0 1^2 0 (2.7) 

\g 0 V 2 J 


( 2 . 8 ) 


For each angle 6 E [0,27r) we define the direction ^{6) := (cos0,sin0). As system (2.1 
is hyperbolic, for each of these directions and a fixed w the matrix 


A(w) = ^iAi(w) + f2A2(w) 


(2.9) 


has real eigenvalues 

Ai = f;-^-c, \2 = v-t, A3 = i7-^+c 

and a full set of linearly independent eigenvectors 



( 2 . 10 ) 


( 2 . 11 ) 


2.3 Lake at Rest 

A trivial, but nevertheless important solution to ( |2.1[ ) is the lake at rest situation, where 
the water is steady and the free surface level is constant, i.e. we have 


z;-(0,0)^andH(f) = Ho. 


From (2.41 we immediately get 


\/h = -Vb 

and therefore (with (2.1 1 - ( |2.3| and F'= (0,0)^) 

0 

0 / .. 


( 2 . 12 ) 

(2.13) 




0 

0 




= -gh 


( 0 

db{x) 

dxi 

db^ 

\ dX2 


(2.14) 


A scheme fulfilling a discrete analogon of (2.14|| exactly is called well-balanced. 
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3 FVEG Schemes 

Finite volume schemes are very popular for solving hyperbolic conservation laws for 
several reasons. They represent the underlying physics in a natural way and can be im¬ 
plemented very efficiently Nevertheless, nearly all of them are based on the solution of 
one-dimensional Riemann problems and therewith a dimensional splitting. This intro¬ 
duces some sort of a bias: Wave propagation aligned with the grid is very well repre¬ 
sented, whereas waves oblique to the grid cannot be caught as accurate. 

In the last decade Lukacova et.al. developed a class of finite volume evolution Galerkin 
schemes, see e.g. [?,?,?]. The FVEG scheme is a predictor-corrector method: In the pre¬ 
dictor step a multidimensional evolution is done, the corrector step is a finite volume 
update. 

In this section we will recall the second order scheme presented in [?]. This method 
will be the starting point for our extensions for computations including dry beds in Sec¬ 
tion Therefore we concentrate on the properties playing a role in this context and limit 
ourselves to the main ideas otherwise. 

3.1 Finite Volume Update 

For our computations, we use Gartesian grids, i.e. we divide our computational domain 
Q in rectangular cells C;, separated by edges E. On the edges, we have quadrature points 
X]^. The subscript i will always refer to a cell, whereas kasa subscript is used as a global 
index for quadrature points. If we talk about the local quadrature points on a single edge, 
we use the index j instead. 

On each cell we define the initial value at as 

u°:^u/(0)«-r^ [ u{x,0)dx (3.1) 

I k; I J Q 

where we use a Gaussian quadrature to approximate the integral. Integrating ( |2.1[ ) on 
each cell, we can then define the update as 

u”+^=u" —T-^ [ ( [ J^{u(x,t))-ndx+ [ S{u{x,t),x)dx\ dt (3.2) 

\Ci\Jt'' \JdCi JCi J 

using the Gauss theorem. Here u” denotes cell average in C/ at time f” and n is the outer 
normal. The solution on the whole domain at time f” is then defined as 

U”(x): = U(f,H)=uf, xeCi. (3.3) 

For an approximation of ( |3.2| , on each edge we define three quadrature points x),; = 
1,2,3, see Fig.j^ These quadrature points are located on the vertices (;=1,3) and the centre 
a = 2) of an edge. The flux over the edge is approximated by using midpoint rule in time 
and Simpson's rule in space, hence we will use the evolution operators from Section [3^ 
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Figure 1: Quadrature points Xj for finite volume update 


to predict point values at the quadrature points at time The flux over an edge is 

then defined as 




n+l/2^ 




■n; 


AflEl 


J^{u{x,t))-ndxdt. 


(3.4) 


At = — t" is the time step, is an approximation to u{xj,t"^^^^) and the ccj rep¬ 

resent the weights of Simpson's rule, i.e. we have ai 3 = g and ^2 = 5 - Finally the source 
term is discretised as 


3 


Si--=g'E‘‘i 


( . y \ 


i-n+l 


Ax At . 


S{u)dx. 


(3.5) 


Here, Ax is the length of the element, hj represents the first component of and 

bj = b{xj). The superscripts stand for the edges surrounding the cell, namely the right, 
left, top and bottom edge. Eqs. (|3.2|-([33]l lead to the fully discrete scheme 


u"+i = u" 


At 

Ax 


E 

E,EcdCi ) 


(3.6) 


The time step is chosen as 

Ax 

Af = umm- 

/■ maxfc|A;t| 


(3.7) 


with A/c the eigenvalues from (2.10 1 and }i<l a CFLnumber. For all the numerical exper¬ 
iments in Section]^ we set /f = 0.5. 


3.2 Evolution Operators 

As mentioned before, we use so-called evolution operators to predict point values of the 
solution for the quadrature points in ( |3.4| |. Indeed, a solution of ( |2.6| can be seen as a su¬ 
perposition of waves. So for a fixed point P = (x,t), we want to identify all the waves that 
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Figure 2: Bicharacteristic decomposition. Left: Bicharacteristic curves for a fixed direction Right: Bicharac¬ 
teristic cone. 


contribute to the solution there. This section will describe shortly the evolution operator, 
present an exact formulation and give an example of a suitable approximation allowing 
an efficient implementation. 

The derivation of evolution operators is based on the quasi-linear form of the system 
( |2.6| . For any given point P, we identify a suitable average value w and linearise the 
system around P: 

Wf+Ai(w)w;,j+A2(w)w;c2 = t- (3-8) 


The waves we are looking for propagate along the characteristics of this system. Thus for 
a fixed direction ^(0), we apply an one-dimensional characteristic decomposition of the 
linearised system. This allows us to identify different wave propagations corresponding 
to the eigenvalues ( 2.10| , the bicharacteristics. The left side of Fig. shows an illustra¬ 
tion. Integrating the decomposed system along the bicharacteristics, we get an integral 
representation of the solution at point P. At this point, the solution still depends on a par¬ 
ticular direction ^(0) and therefore does not respect waves coming from other directions. 
Thus we perform the decomposition for all angles 0 G [0,27r) and average the solution at 
P over 0. This yields the exact evolution operator of ( |3.8| . The combination of all bichar¬ 
acteristics yields the bicharacteristic cone shown in the right picture of Fig. We introduce 
the following notation for the peak P = (T,f” -|-t) and points on the sonic cone: 


Qo:={xi + TV-i,X 2 + TV 2 ,t") (3.9) 

Qo:={xi + {t" + T — t)vi cos 9 ,X 2 + {t’^ + T—t)v 2 sin0, t) (3.10) 

Q:= (vi-|-T(c-|-tii)cos0,V2-|-T(c-|-ti2)sin0,f") (3.11) 

Q:= {xi + {t" + T — t){c+Vi)cos 6 , X 2 + {t" + T — t){c+V 2 )sm 6 ,t) (3.12) 


Qo is the centre of the sonic circle at time t = f”, Qo denotes a point on the inner bicharac¬ 
teristic connecting P and Qo, Q is a point on the perimeter of the sonic circle at time t = t" 
and Q denotes a point on the mantle of the sonic cone at an arbitrary time Fg [t",t"+T]. 
After some tedious calculations, see e.g. [?], we get the evolution operators for the 
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SWE: 




c2n 


0 




2n, 


HQ) - -ivi{Q)cose+v 2 {Q)sme) de 
s 

ff+T /'In 

/ {bx-^{Q)cos9 + bx2{Q)sme)dedt 


(3.13) 


I't+T 


c 


t + T-tJo g 

(•'In 


-(^i(Q)cos0+Z72(Q)sin6) dOdt 


1 1 O' 

Vi{P) = -Vi{Qo) + — / -3^(Q)cos0+z;i(Q)cos^0+z;2(Q)sm0cos0d0 

z In Jo c 


rt+T 


-i 

2nJt 

1 


(Qo) + ^xi (Qo) dt 

f+T /'In 

/ (&xi (Q) COS^ 0 + &X 2 (Q) sin 0 COS 0) cf 0£f f 

20 

J />27r 


(3.14) 


1 + T — t Jo 
1 


t^ 2 (E) = ^t^ 2 (Qo) . „ , 

2 In Jo 


(i2i(Q)cos20+?72(Q)sin20) (i0(il 
— ^? 2 (Q)sm 0 +z;i(Q)sm 0 cos 0 +r’ 2 (Q)sin^ 0 d 0 


_g 

2 if 

271. 


rf+T 


(Qo) +^X 2 (Qo) 

ff+T />27r 


(3.15) 




(fc^i (Q)cos0sin0 + &;,;2(Q)sin"^0) dOdt 

1 

(771(0)81020 + 172(0)00820) d 6 dt. 


1 + T — t Jo 


For efficient computations these operators have to be simplified. In [?], the authors follow 
[?] and present approximations of (3.131 - ( 3.15| that provide exact solutions of some one¬ 
dimensional Riemann problems. For piecewise constant data, these approximations read 




pin 




H(Q)--(i7 i(Q)sgn(cos0)+i72(Q)sgn(sin0)) 

T f'2.7Z 

{vibx,{Q)+V2bx2{Q))d9, 


d 9 


^^(^) = ^Io [-|w(Q)sgn(cos0) 




pin r 


+17i(Q) ( cos^0 + - ) +772(Q)sin0cos0 


-|H(Q)sgn(sin0) 
+77i(Q)sin0cos0+772(Q) ( sin ^0 + ^ 


d 9 , 


d 9 . 


(3.16) 


(3.17) 


(3.18) 
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Figure 3: Left: Intersection of the sonic cone at quadrature points with grid cells. Right: Stencil of a quadrature 
point 


The corresponding operators for piecewise (bi-)linear data are given as 
h{P) = H{Qo)il-^)-b{P) + y^" H{Q)de 

C 

~JnJo (^i(Q)‘^°®^+^2(Q)sin0)d0 
+ {vibrAQ)+^2bx2{Q))de, 

v^{P) = viiQo){l-^) + ^ rH{Q)cosede 

4 CTT JO 

1 

+ iy [!7i(Q)(l+3cos^0)+3z72(Q)srn0cos0] dO, 
Vi{P)=V2{Qq){ 1-"^) + ^ f^"HiQ)smede 

4 CTT Jo 

1 

T-y [3!7i(Q)sin0cos0+!72(Q)(l+3sin^0)] 


(3.19) 


(3.20) 


(3.21) 


We introduce the operator E(W) as a shorthand for evaluating these evolution operators 


at all quadrature points Xk for any given numerical data W defined analogously to (3.3'. 
We will refer to the operators for piecewise constant data ( 3.16|l- ('3.18|l as E“'^*(P) and 
Ebiim(p) -will denote the operators for piecewise bilinear data (3.19 - (3.21 1 . 

In our schemes, these operators are evaluated at the quadrature points vj- of the finite 
volume update defined in ( |3.4| . Thus all data contributing to the evolved values is de¬ 
rived from the cell values next to the quadrature point. We therefore define the stencil S;t 
of a quadrature point as 

Su:={Ci\xkedCi}. (3.22) 

An example of the intersection of the cone with grid cells and the resulting stencil is 
shown in Fig.|^ The suitable average value wjt used in (3.8 1 is chosen as 




|S^ 


LQeSj: 


(3.23) 
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We also tried a local Lax-Friedrichs update at the prediction points to get a better lineari¬ 
sation. The numerical results were almost exactly the same, so we chose the averaging 
procedure (3.23| for our computations. 


3.3 Numerical Representation of the Bottom Topography 

For finite volume schemes, the numerical representation of the bottom topography plays 
a crucial role in well-balancing as well as positivity of the scheme. In [?] Audusse et.al. 
use cell averages of the bottom for the computation of the free surface and reconstruct 
the free surface and the water height. The reconstruction of the bottom then results as 
the difference between slopes of H and h. In [?] Kurganov and Petrova propose to use 
a piecewise linear approximation of b instead of b itself by taking the values of b at cell 
corners. The cell average of b is then computed as the average of the corner values. 

For the FVEG schemes it is necessary to define some value of b not only for cell aver¬ 
ages and the reconstructed slopes, but also at the quadrature points where the evolution 
operators are evaluated. There is some freedom in doing this, as the source term dis¬ 
cretisation ( |3.5| respects the well-balancing property independently of the reconstructed 
slopes of the bottom topography. As the evolution operators for the water height com¬ 
pute the free surface first and derive the actual water height via h{P) = H{P) — b{P), the 
only necessary condition for b{P) is b{P) <H{P). 

In this work, we will define the cell averages of b as in 0- For the quadrature points 
on cell corners, we set 

— E bi^b{xk). (3.24) 

I i-CiESit 

The values of at the centres of each edge are linearly interpolated from the neighbour¬ 
ing corners. While the latter condition has been derived in [?] to ensure well-balancing 
on adaptive grids, the formula for the comer points will turn out to be helpful for the dry 
bed case. 


3.4 A Multidimensional Entropy Fix for the FVEG scheme 

It is well known that the weak solution of a Riemann problem for conservation laws is not 
always unique, and an entropy condition is needed to single out the physically correct 
solution. This has its correspondence on the discrete level, where conservative numerical 
schemes may converge to entropy-violating solutions. This notorious difficulty seems 
to appear only near sonic rarefaction waves, where the flow changes from subcritical to 
supercritical velocity [?]. Various researchers have proposed so-called "entropy-fixes" 
for numerical schemes. In particular, we would like to mention Harten's and Hyman's 
entropy fix for the the Roe solver [?] (see also the discussion in [?]). 

The FVEG schemes considered here make no exception, and may compute entropy 
violating solutions, see the left picture in Fig. As is well known for classical finite 
volume methods, this effect is less visible (though still there) for second order schemes [?]. 
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Figure 4: ID dam break problem, solved with first order FVEG method. Left: solution without entropy fix. 
Right: solution with entropy fix from [?] 




Figure 5: Position of sonic cones for the discrete sonic rarefaction with subsonic U; and supersonic u,. Left: 
cone for 0 / < C;. Right: cone for w, > c. 


In order to make our point clear, we therefore focus on first order computations for the 
rest of this section. 

In [?], Lukacova and Tadmor proposed an entropy conservative variant for rarefaction 
waves computed by certain Riemann solvers, see also [?]. They applied this technique 
successfully to the finite volume corrector step of the FVEG scheme. They derived just the 
right amount of viscosity that one should add to the scheme to fulfil the entropy equality. 
Fig. [^clearly shows the effectiveness of the scheme: While the standard FVEG scheme 
produces an entropy violating shock, the entropy conservative scheme clearly reproduces 
the correct rarefaction wave. Nevertheless, the scheme from [?] does not appear perfectly 
suitable for our needs. First, the proposed fix requires the characteristic decomposition of 
the jump of the conserved quantities across an edge. As the decomposition is not needed 
for the FVEG schemes, this is an undesired computational extra cost. In the context of 
dry boundaries, we should also mention that the decomposition matrix becomes very ill 
conditioned when /j—f 0. The second point is that the scheme from [?] has been developed 
for the one-dimensional case. Although it can be applied dimension wise, this approach 
somewhat spoils the multidimensional spirit of the FVEG methods. 

We therefore propose a new approach to solve the entropy problem. It is not based on 
a flux correction, but on the correct evaluation of the EG operators. To motivate our solu¬ 
tion, we take a closer look on a discrete one-dimensional Riemann problem that should 
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Figure 6: Position of sonic cones for the discrete sonic rarefaction with subsonic U; and supersonic u,. Left: 
cone for averaged value u. Right: superscribed cone. 


result in a transonic rarefaction, i.e. the flow is subsonic in upwind direction and super¬ 
sonic in downwind direction. Thus let us assume we have two adjacent cells C/ and C,- 
with cell averaged data 

u/ = {hi,vi,0), Vi < Cl and u,- = {hr,Vr,0), v,- > Cy. (3.25) 


Here c is the speed of sound defined in (2.51. To evaluate the evolution operators, we 
start with the sonic cones defined in ( |3.9| - ( |3.12 |. For simplicity, we limit ourselves to the 
quadrature point at the centre of the edge. The sonic cones resulting from u/ and are 
sketched in Fig.|^ We can see that the cone resulting from is located completely in C/. 
Depending on the exact values of u/this can also be the case for the sonic cone result¬ 
ing from the averaging procedure ( 3.23| |. In other words: We use an evolution operator 
resulting from a supersonic linearisation in a regime that is subsonic. At least for the first 
order operators (3.161 - ( 3.18| , this means that the predictor step exactly reproduces u/, 
which is then used for the flux evaluation. This corresponds to the generalised upwind 
method which is known to compute entropy violating solutions in some situations, cf. 
e.g. [?]. 

Thus the core of the problems seems to be the wrong domain of dependence for the 
predictor step. In case of a sonic rarefaction, the sonic cone should always include both 
regions, the subsonic as well as the supersonic one. As this is not guaranteed, we modify 
our method by extending the sonic cone if necessary. In a transonic situation we drop the 
sonic circle resulting from the averaging procedure ( 3.23| . Instead, we use a circle which 
comprehends all the circles defined by the cell averaged values in the corresponding 
stencils, see Fig. for an illustration. The exact formulation used for our schemes is as 
follows. Given two circles with midpoints v, and radii r„ i — 1,2, we compute the new 
circle as 


ri+r2+d 

2 

x = xi + {r-ri) ^'^ 

with d = 11X 2 ~ H11 2 - However, if one circle comprehends the other one, we just choose the 
bigger circle as our new one. If the stencil of the evolution point consists of four cells, we 
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Figure 7: ID dam break problem, solved with first order FVEG method. Left: solution with entropy fix from [?] 
Right: solution with new entropy fix 


apply the same formula for the neighbours on the two diagonals first and then again for 
the resulting circles. 

In Fig. 1^ we compare the results of the entropy stable scheme from [?] and our new 
approach. They both clearly solve the entropy problem, with a slight advantage of the 
scheme by Lukacova and Tadmor. However, as we previously pointed out, the new 
method is more efficient. Compared to the scheme without any entropy fix, the new 
one took only about 2% extra time, whereas the approach from [?] needs about 8% more 
computational time. As the new scheme is also suitable for computations including dry 
areas, we chose it for all the numerical experiments in Section]^ 


4 Dry Bed Modifications 


To extend our schemes to computations including dry beds, we have to guarantee two 
properties: the positivity of the water height, and the well-balancing under the presence 
of dry areas. In literature, this is mainly achieved by two basic ingredients: a positivity 
preserving reconstruction and an additional time step constraint. Examples can be found 
in [?,?,?]. 


For the FVEG schemes, these measures fall short of the aims. The additional predic¬ 
tor step via the evolution operator prevents a direct proof of a positivity property. One 
reason for this is the extended stencil: the flux over an edge is computed using more cells 
than the direct neighbours (see Fig.[^and|^. Another problem is the complex evaluation 
of the operators and with them the flux which makes an analysis of the positivity at least 
challenging if not impossible. 

Regarding the well-balancing, a sophisticated reconstruction is not enough. From 
(|3.5| it is obvious that the reconstruction does not directly affect the balancing of flux and 
source terms. The core of the problem is that the lake at rest described in (|2.12 1 changes 
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to 

z/-(0,0)^andH(x)-|^° Mx)>0 

lo(xJ else 

if dry areas are included. Thus for our schemes, we have to find evolution values for the 
water and bottom height that can handle properly the occurrence of this situation, i.e. 
which avoid the generation of spurious waves at the shoreline. 

In this section, we will present an alternative approach to ensure the positivity of 
the water height as well as modifications of the finite volume update and the evolution 
operator to ensure the well-balancing property. We make sure that the changes do not 
affect the scheme away from dry regions. 


4.1 A General Positivity Preserving FV Update 

For the derivation of a positivity preserving scheme, we study the first component of the 


finite volume update (3.6 1 


At 




E,EcdCi 


where 




(4.2) 

(4.3) 


is the first component of the flux vector the following, we will modify the flux 

in order to guarantee positive water height. The technique presented here is applicable 


to arbitrary finite volume fluxes, and the FVEG flux given in (3.41 below is only a special 
case within this framework. 

The basic idea of our method is to cut off the outgoing fluxes as soon as all water 
which has been contained in the cell at the beginning of the time step has left the cell via 
the outgoing fluxes. We will call this time the draining time. For convenience, we will 
later rewrite this as a reduced time step Af£ on these edges, but in fact the finite volume 
update will always advance the solution by one global time step At. 

Thus our first step is to separate the fluxes contributing to the outflow off a cell C 
from those contributing to the inflow by setting 


:=max{7^E,0} 

:=min{'KE,0}. 


This allows us to rewrite the update (|4.2| as 


E,EcdCi E,EcdCi 




outflow 


inflow 


(4.4) 

(4.5) 


(4.6) 
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Now we introduce the draining time by 


Atr. 


C;,drain • 


Axhf 


^E,EcdCi'^E 


(4.7) 


Once again, at time Ate. drain all water which was originally contained in cell C, has 
flown out, so 


hf 


AIq, drain 

Ax 


E «E=o- 

E,Ec3Q 


(4.8) 


outflow 


Suppose now that Atq drain < At. For t G [t”+Atq drain/no more water can leave the 
cell, at least not water which was originally contained in cell C,. Therefore we assume that 
there is no outgoing flux for times beyond the draining time, and introduce the cut-off 

flux by 


nUt)-= 



for f" < f < t"+Afc, 4 rain 
for t ^ t TAIq drain 


(4.9) 


Now we integrate the draining flux in time and obtain the cut-off (or draining) finite 
volume flux 




r+ 


'HE{t)dt = 


rfl'+Atc drain 


It” 


7^+(u”)dt = At, 


Ci.drain 


ntivL^ 


(4.10) 


Before introducing the positivity-preserving modification of the finite volume scheme 
( |3.6| l, we rewrite the cut-off in the flux as a local cut-off in the time step. This cut-off time 
step is defined for each edge E and takes into account the upwind cell (E): 


AtE :=min ( At, At 


(E),drain 


(4.11) 


Now we replace the finite volume flux 7E in ( |3.6[ ) by the cut-off finite volume flux defined 
in ( 4.10 1. Using ( 4.11| |, this leads to the modified general update (3.6 1 




Ax 


(4.12) 


E,EcdCi 


Theorem 4.1. Assume we have a conservative finite volume scheme for the solution of the shallow 
water equations that can be written in the form \3.6\. Then the modified finite volume scheme 
( 4.12[ ) with locally cut-off flux ( |4.10| |) (respectively locally cut-off time step ( |4.11[ )) is positivity 
preserving. 

Proof. From the first component of ( |4.12[ ), combined with definitions ( |4.11| | of the cut-off 
time step AtE and ( |4.7| | of the draining time Afc„dram/ the water height at the new time 
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step can be bounded as follows: 




E,EcdCi 


Az 


>K- L 

E,E(ZdCi ^ 


> K - £ 

E,EcdC 


At, 


C (E),drain 

Ax 


m 


hf 


= 0 


drain 


E 


L-l '^E 

E,EcdCi 


□ 


Remark 4.1. The time step Afg used in the new finite volume update ( |4.12 | and defined 
in ( 4.11| | might seem to be local for each edge. We would like to stress, however, that the 
finite volume scheme ( 4.12| | still advances the solution by one and the same global time 
step Af. The apparent contradiction is resolved by considering equation ( 4.10| |: the time- 
integral of the flux is still over the global interval However, the flux "K 

cut-off in the presence of vacuum, see 


^E,drain 


IS 


4.2 Well-balancing at the Shoreline: the Finite Volume Update 

In the derivation of the positivity preserving finite volume update, we so far neglected 


the source term. Its introduction to the new scheme (4.121 rises the question which time 


step should be used for the source term. To maintain the well-balancing, the source term 
and the gravity driven parts of the flux must be multiplied with the same time step. This 
is in contradiction to the definition of Af£, which may change for different edges of the 
same cell. On the other hand, the reduced time step is not necessary for the momentum 
equations. We therefore shift the gravity driven components of JF into the source term, 
i.e. we define 





hV2 \ 

hviV 2 I and <S*(u,T) 


hvl 



By replacing JF with F"* in (3.41 and changing (3.51 to 




/=i 


0 






(4.13) 


(4.14) 
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we can rewrite (3.6 1 as 


u"+i=uf-^ 
' ' Ax 


AtEJ^*^+AtS* 


KE,EcdC, 


(4.15) 


This formulation ensures the well-balancing of the scheme even in the case of modified 
local time steps. Away from the shoreline we have AtE = Af V£ and ( 4.15| l equals the 
original update (|3.6|. 


Remark 4.2. The rearrangement of the advective and gravity driven parts of the equa¬ 
tions in (4.13| is independent of the scheme. Thus in principle, every finite volume 
scheme can be reformulated as in (4.151. To obtain a well-balanced scheme, we only 
need a discretisation of S* analogously to (4.141 that preserves the lake at rest solution in 
the presence of dry boundaries. 
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H 


B 


H 


B 


Figure 8: Lake at rest with dry boundaries. Solid line: free surface, dashed line: bottom topography, filled 
circles: H at evolution points, empty circles: B at evolution points. Left: Real situation. Right: Numerical 
Representation with evolution values. 


4.3 Well-balancing at the Shoreline: the Evolution Operator 

The lake at rest situation with dry beds described in (|4.1| is only preserved if the numer¬ 
ical flux and source terms in ( 4.14| are exactly balanced, or equivalently if the evolution 
operators reproduce the lake at rest. This is not necessarily the case if the stencil of a 
quadrature point contains dry cells, as is demonstrated in Fig.|^ If bi> Ho for a cell in the 
stencil, the resulting bottom value from (3.24|| can be higher than the free surface in the 
wet cells. Using the approximate evolution operators ( 3.16| and (3.191, it is easy to see 
that in the lake at rest case the evolved water height is also positive, leading to an even 
higher free surface at the quadrature point. In this case the combined flux and source 
term S* from ( 4.14[ | does not vanish anymore and introduces unphysical waves starting 
from the dry boundary. 

To avoid the creation of these waves, we modify the data used for the predictor step 
at the interface. First, we replace the stencil Sjt by S| defined as 


Sl:={Ci\CieSkAhi>0} 

which allows us to determine the maximal free surface level at xjt as 

Hjt = max(H,). 

C* 


Now in ( 3.23| and for the evaluation of the evolution operators we set 

if Q ^SlAbi> Hu. 


(4.16) 


(4.17) 


(4.18) 


In all other cases, we leave the values unchanged. The modification, which is illustrated 
in Fig. en sures that the free surface is correctly represented in the source term compu¬ 
tation ( 4.14| . We also avoid an unphysical flooding of mounting slopes. In [?, ?] a similar 
technique is used on triangles. 

Finally, even in the presence of wet, but nearly dry cells in the stencil, the expressions 
for h{P) in (3.16' and (3.19| can become negative if hi is small in the surrounding cells. 
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H 


o 


B 


o 


X 


Figure 9: Lake 
like in Fig. 


at rest with dry boundaries, modification |4.18| for computation of evolution values. 


Symbols 


This cannot necessarily be cured by a smaller time step, as substantial parts of the ex¬ 
pressions do not depend on At. With this restriction in mind, we propose the simplest 
solution: Whenever we have h{P) <0, we set h{P) =Vi{P) — V 2 {P) =0. 


4.4 Treatment of Nearly Dry Cells 

In our schemes, we consider a cell Ci to be dry when hi < Eh, where we have chosen 
Eh = 10^^. In dry cells we set hj = w, = Vi = 0. A well known problem occurs when hj is 
close to that value: the velocity v = hv/h can become singular due to small numerical 
errors in the conserved variables. This leads to very small time steps which in the end 
can basically stop the computation. This problem has been discussed e.g. in [?,?], where 
different strategies have been proposed. In [?], Kurganov and Petrova propose to de- 
singularise v by multiplying it with a certain factor / < 1 whenever h falls below a certain 
threshold Ev In [?], the authors just set z? = 0 whenever h < Ev 

In this work, we use a different approach. As the solution at the dry boundary is 
always a rarefaction wave, the flow velocity will grow smoothly when water floods for¬ 
merly dry areas. We will therefore limit the velocity in nearly dry regions depending on 
the velocity in flooded areas. We define the reference speed 

2 (4.19) 

Xi-Xj\\oo. (4.20) 

Whenever we have hi < Ev and ||zf;|l > v^gf, we set the new velocity to 

V*=Vref(2-^y (4.21) 

such that ||zi|| is smoothly limited to a value between v^ef and 2Vrgf. The velocity compo¬ 
nents in Ci are then defined as 

Vi = v*d (4.22) 


where 


Vygf\= max ||!7| 

V.hi>£y 


Ax 


-•ref 


, Lygf :=max 


h] 
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with d the unit vector pointing in the same direction as the vector of discharge {hv\,hv 2 Y■ 
This approach appears us to be a better representation of the physics of the flow, as the 
velocity at the front is not necessarily vanishing. 


4.5 The FVEG Algorithm 

Before summarising the whole FVEG algorithm, we will spend a few words on the re¬ 
construction needed to evaluate the evolution operators for piecewise linear data ( 3.19| - 
(3.211. As the operators are computed from the primitive variables w, these are a natural 
choice for the reconstruction Thus in each cell we need the linear function 


RAx (W,) (x ) IQ : ^ w/ -6 V W/ • (f - f;) -h (w;) (^- V') 1 (^- V-) 2 • 


(4.23) 


The derivatives (w;)xj,(wj);i ;2 {vii)xiX 2 ^re computed from the slopes between cell av¬ 
erages, cf. [?]. In this paper, we use the continuous, piecewise bilinear recovery described 
in [?] without any limiters. The piecewise bilinear functions are uniquely defined by the 
averages at the cell corners, which are already computed for the evaluation of the evolu¬ 
tion operators, see ( 3.23| . Then the w from (4.231 is exactly the average of the averages at 
the cell corners, thus the resulting reconstruction is not necessarily conservative. We will 
therefore use the combined evolution operator 


E(W): = (Rax(W) ) -6E™’'®*(W- W), 


(4.24) 


where the first order correction E“”®*(W —W) is necessary for stability, see [?] for an in- 
depth discussion of this issue. 

Although the conservative correction introduces some oscillations at steep fronts, the 
piecewise bilinear reconstruction has a clear advantage in the vicinity of dry areas. As 
the averaged water height at the cell corners is non-negative via eqs. (4.161 - ( 4.18| , the 
reconstruction is also non-negative by design. In dry cells, we set 


= 0if hi = {). 

We refer to [?,?] for further details concerning the reconstruction strategy. 
The complete algorithm now reads as follows: 


(4.25) 


Algorithm 4.1. 1. From given conservative data u" and hi at time f”, compute the 

nonconservative variables H",v" ^,^2 

2 . apply the reconstruction operator Rax to Hf,v\ , arid 

3. compute the evolution operators 

4. evaluate the advection fluxes from (4.13|l 


5. compute the gravity driven flux and source terms S* from (4.14 
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Figure 10: Circular dam break over dry bed, solution at t=1.75s. Left: 3D view. Right: 30 contour lines 
between 11 = 10.2 and 11 = 0.2 



6 . perform the finite volume update ( 4.15| l 

We will finish the section by a proof of the well-balancing property of the scheme. 


Theorem 4.2. Suppose that we have a numerical solution respecting the lake at rest solution {4.1^ 
with dry boundaries. Then the FVEG scheme { 4.15} together with the modifications described in 
Section [4J]preserpes this state. 


Proof. For the lake at rest state with v= (0,0)^, the advective parts of the fluxes IF defined 
in (|2.2| and IF* defined in (4.13' are all zero. Thus for all edges we have Af£ = At and the 
original finite volume update ( 3.6| and the modified update ( 4.15| are the same. Then 
Theorem 2.1 from [?] states that the scheme is well-balanced provided that the predicted 
point values used for the flux evaluation also satisfy the lake at rest situation. 

We will now show that all data used for the reconstruction and for the evaluation 
of the predictor step satisfies the requirements of Theorem 3.1 from [?]. From definitions 
(4.161 and ( 4.17| we see that for all evolution points the averaged free surface is computed 
as Hq, which is the free surface level in all flooded cells. As all velocities are zero by (4.11 
respectively ( |4.18 l, the averaging also returns zero. Finally, the reconstruction procedure 
is based on the averaged point values and therefore in all flooded cells we have w^j = 
Wx 2 = = 0- Iri dry cells we have the same result by definition ( 4.25| . Thus we can 

apply Theorem 3.1 from [?] and this concludes our proof. □ 


5 Numerical Results 

5.1 Dam Break over Dry Bed 

This is a classical test case where we simulate the complete break of a circular dam sep¬ 
arating a basin filled with water from a dry area. The computational domain is [0,100]^ 
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Figure 11: Circular dam break over dry bed, solution at t=1.75s with r=||x|| 2 . Left: free surface. Right: 
velocity 


and we set Az = 1, the water filled basin is located at r = ||x|| < 60. In the basin we set 
Ho = 10 and elsewhere Hq = 0 and the initial velocity is vq = (0,0)^'' in the whole domain. 
Reference solutions can be found in [?,?,?]. 


In Fig. 10 we see a 3D view and contour lines of the water height. Fig. 11 shows the 


water height and velocity at different lines through the domain. The resulting rarefaction 
wave is almost perfectly symmetric, the oscillation due to the reconstruction strategy 
is restricted to three percent of the water height (we have maXjHj = 10.269). We see a 
small bump at the drying wetting front, but the front position and velocities are well 
represented. Thanks to the new entropy fix, there is no unphysical shock visible in the 
transcritical region. 


5.2 Wetting/Drying on a Sloping Shore 

This test case was proposed by Synolakis in [?] and computed in e.g. [?,?]. It describes 
the run-up and reflection of a wave on a mounting slope, with the initial solution given 
as 


Ho{x) =max(/,&(T)), i;o(^) = ^Ho{x),0^ 


where 


f{x) = D+dsech^{'y{x-i-Xa)). 

As in [?, ?], we set D = l,ci = 0.019 and 

T=Y|^. Xa = ^^arcosh{V20). 
For the bottom topography we have 

0 < 2Xa 

b{x)=b{xi)=i Xi-2Xa 


19.85 

22 
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The computational domain is Q = [0,80] x [0,2] and the grid size Ax = 0.04. We prescribe 
open boundary conditions in Xi direction and periodic ones in Xz direction. 

In Fig. we present the water height during the run-up and drying process together 
with the analytical solution. Details on how to obtain the analytical solution can be found 
in [?, Section 3.5.2]. At time t = 9, the wave has almost reached the shoreline, and shortly 
after, at f = 17, the wave reaches the maximal run-up on the shore and the drying pro¬ 
cess starts. During the run-up, the agreement with the analytical solution is excellent. 
The drying process is reproduced very accurate, too, although compared to the run-up 
process, there are small deviations between numerical and analytical solution with re¬ 
spect to the minimum of the water height at f = 23. At t = 28, where the reflected wave 
starts leaving the domain, these deviations persist, but stay very small. During the whole 
simulation, no oscillations or other perturbations at the wetting/drying front are visible, 
demonstrating the well-balancing capabilities of the scheme. The scheme also returns 
quickly to the lake at rest solution presented in the last picture (f = 80). where the water 
surface is almost flat. All in all, the results are very satisfying and compare well to the 
other numerical solutions presented in [?,?]. 


5.3 Vacuum occurrence by a double rarefaction wave over a step 

To test how the scheme handles the drying of formerly flooded areas, we compute a test 
case proposed in [?]. It describes two separating waves over a non-flat bottom. The 
computational domain is a pseudo-ID channel given as Q = [0,25] x [0,0.5]. We set the 
initial free surface height to = 10 and the discharge and bottom topography to 


hv-i{x,0) 


350 ifxi>50/3 

—350 else 



if25/3<xi<25l2 

else 


(5.5) 


Like in [?], the computation is performed on a grid with 300 grid cells in Xi direction. 

In Fig. 13 we show the free surface and the discharge of the solution at different times 
t. At time t = 0.5, several waves have arisen from the interaction of the supersonic flow 
with the bottom topography. Due to the reconstruction strategy described in Sec. 4.5 the 
rarefaction waves are smoothed out at the bottom and show small oscillations at the top. 
However, for the following times we see an accurate representation of the waves flowing 
over the hump (f = 0.25) and leaving the domain (f = 0.45,0.65). No spurious modes are 
introduced during the drying process, neither in the flat region nor at the edges of the 
hump. 


5.4 Thacker's Periodic Solutions 


We present two exact solutions of (2.11 proposed by Thacker in [?]. They both describe 
oscillations of a free surface in a parabolic basin with a free shoreline. The basin is defined 
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Figure 13: Vacuum occurrence, solution at different times. Left: free surface. Right: discharge. 




Figure 14: Thackers curved solutions. Left: t = T. Right: t = 3T 


as 


b{x) = b{rc) = -Ho 



(5.6) 


Tc defines the distance from the basin's centre, Hq the height of the centre and a is a 
parameter. We will define two functions that describe solutions of (2.1 1 with 

h{x,t)=max{f{x,t)—b{x),0). Both test cases will be computed on the domain Q=[—2,2]^. 
The presented results have been performed with a grid size of Ax = 0.4. Reference solu¬ 
tions can be found in [?, ?, ?]. 


5.4.1 Thacker's Curved Solution 

The first function results in a curved oscillation over b, it reads 


f{rc,t) = Ho 



Vl-A^ 

1 —Acos(a;f) 



A. 

{1 — Acos{cot))^ J j 


(5.7) 
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# cells 

L“ 

EOC 


EOC 


EOC 

25 X 25 

9.8038e-03 


1.4526e-02 


9.2884e-03 


50 x50 

3.6191e-03 

1.44 

3.9038e-03 

1.90 

2.7762e-03 

1.74 

100 xlOO 

1.5252e-03 

1.25 

1.3127e-03 

1.57 

9.5937e-04 

1.53 

200 x200 

1.1820e-03 

0.37 

4.6649e-04 

1.49 

3.8549e-04 

1.32 

400 x400 

5.3221e-04 

1.15 

1.7806e-04 

1.39 

1.4907e-04 

1.37 


Table 1: Experimental order of convergence (EOC) for Thackers curved solution. Error in water height in 
different norms. 


Fceiis EOC ? EOC I? E^ 


25 X 25 

3.3855e-02 


4.6507e-02 


2.7148e-02 


50 x50 

1.7455e-02 

0.96 

1.8179e-02 

1.36 

1.1660e-02 

1.22 

100 XlOO 

1.0543e-02 

0.73 

1.0486e-02 

0.79 

6.6938e-03 

0.80 

200 x200 

8.2376e-03 

0.36 

8.3640e-03 

0.33 

5.4658e-03 

0.29 

400 x400 

7.1238e-03 

0.21 

7.8559e-03 

0.09 

5.1747e-03 

0.08 


Table 2: Experimental order of convergence (EOC) for Thackers planar solution. Error in water height in 
different norms. 


Here, co = VspoA? is the frequency and for a given ro>0, A is the shape parameter 


A = 


2 2 
« -^0 

a^+r, 


0 

Eor the computation we set a = 1, Hq = 0.1 and tq = 0.8, which results in an oscillating 
period of T ^2.22. The initial velocity is set to Vq = (0,0)^'". 

5.4.2 Thacker's Planar Solution 

The second solution is a planar surface rotating around the basin. The corresponding 
function is 

(5.8) 


= +2{x-xc) ■ {cos{(vt),sm{(vt)y’') 


with cv = a?- the frequency and rj another parameter. Here, we set a = 1, Hq = 0.1 

and 7/=0.5. The resulting period is then T «4.44. The initial velocity in the wetted domain 
is given as vq = (0,7/a;)*''. 

We present the water height after one and three oscillations along the line X 2 = 0 in 
Eig. 14 for the curved solution and in Eig.j^for the planar solution. The exact solution is 
very well reproduced, independent of the shape of the initial solution. We can see a slight 
smearing after three periods for the curved solution, where the maximum value at the 
centre is reduced and the drying/wetting interface has been pushed outward. Similarly, 
the interface of the planar solution has also moved a little bit inwards at f = 3T. Again, 
for both solutions there is no production of spurious waves at the dry boundary. 
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Figure 15: Thackers planar solutions. Left: t = T. Right: t = 3T 



Figure 16: Circular island, lake at rest situation at t = 5. Left: whole domain. Right: zoom on island. 


In Tables]^ and 1^ we present a convergence study for the two test cases. The experi¬ 
mental order of convergence for the curved solution is well better than one, which meets 
our expectations. The errors are slightly better than in [?]. For the planar solution, how¬ 
ever, the order quickly drops to zero. The problem seems to raise from the boundary of 
the wetted domain, where we have supersonic velocities tangential to the bottom slope. 
So the problem might be related to the evolution operators, as they produce inexact solu¬ 
tions in other supersonic situations as well, see the discussion in Section 3.4 and|^ 


5.5 Wave Run-up on a Conical Island 

In this case we simulate the run-up of a solitary wave over a conical island. It has been 
performed experimentally at the U.S. Army Engineer Waterways Experiment Station, 
see [?,?]. The computational domain is Q= [0,25] x [0,30] and we set Ax = 0.2. The centre 


27 























January 16, 2015 






eu 

4.44089e-16 

6.50361e-14 

3.14147e-15 


2.69215e-15 

2.53222e-13 

1.24917e-14 


2.92617e-15 

2.73056e-13 

1.33148e-14 


Table 3: Lake at rest around a conical island. Errors at t = 5. 



Figure 17: Run-up on a circular island. Left: wave approaching island, t = 7.9. Right: run-up at front of the 
island, f = 9.1. 


of the island is located at Xq = (12.5,15) and with r = ||x — fc || its shape is given by 


b{r) = < 


0.625 

(3.6-r)/4 

0 


r<l.l 

r<3.6 

else. 


(5.9) 


The initial free surface is given by Hq = 0.32. We start by giving an example of the well¬ 
balancing capabilities of the scheme and compute the lake at rest situation until t=5. The 
results are shown in Fig. The lake at rest is perfectly preserved, which is confirmed 
by the errors given in Table For all the 3D views of this example, the vertical axis 
representing the free surface was scaled by a factor of 25 to emphasise the results. 

We now compute the actual wave where at time f=0 a wave enters the computational 
domain at Xi = 0. The height of the wave is given by 


= Ho+ciHo 


\cosh{^y/gHo/L{t-3.5) 


with L = 15, a = 0.1 and ^ \ 3ci{l+a)L'^/ (4Hq), cf. [?,?]. We present 3D views of the 


solution in Fig. 17 18 and 0^ Fig. 17 shows the wave approaching the island and the 


instant of its maximal run-up. Here, as in all the following figures, the region in front 
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Figure 18: Run-up on a circular island. Left: lateral run-up, f = 10.7. Right: symmetric waves around the island, 
t = 12.1. 



Figure 19: Run-up on a circular island. Left: run-up behind the island, t = 13.7. Right: reestablished lake at 
rest at f = 40. 


of the wave remains completely unaffected, once again confirming the well-balancing of 
the scheme. Due to the run-up, the wave is slowed down at the island, finally resulting 
in a reflection of the wave presented in Fig. This leads to the formation of two sym¬ 
metric waves surrounding the island with a noticeable run-up also on the sides of the 
island. Away from the island, we observe the circular reflected wave approaching the 
boundaries of the computational domain. Behind the island, the lateral waves reunite 
and produce a second peak in the run-up on the lee side of the island, see Fig. 19 Finally, 
the wave leaves the domain and the surface returns to the lake at rest, as shown in the 
last picture. 


Moreover we show the evolution of the free surface at chosen points in Fig. The 
position of the gages is given by = (6.36,14.25), X(, = (8.9,15), xg = (9.9,15), Xig = 
(12.5,12.42) and X 22 = (15.1,15). For comparison we also display the measured data ob- 
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tained from [?] as well as the results from the residual distribution (RD) scheme presented 
in [?], which have been computed on a comparable gricQ Both numerical schemes pro¬ 
duce steeper wave fronts than the experiment, which results from the lack of dissipation 
in the shallow water model. The FVEG scheme shows a stronger steepening than the 
RD scheme, and also predicts slightly higher waves at all gages. At gages 3 and 6, the 
run-down process is pretty similar with both schemes, while the FVEG scheme returns 
quicker to a constant water height. Both schemes produce less accurate results at gage 9. 
The FVEG line is shifted to the upper left, thus resulting in a less pronounced run-down 
and again it returns quickly to a constant height. The RD scheme provides smoother re¬ 
sults with a better approximation of the maximal run-down, but stays below the original 
water height for a longer period. At gage 16, the FVEG scheme gives a quite accurate 
representation of the wave, where the RD scheme introduces an undershoot after the 
first wave. At the last gage, the graph of the FVEG scheme is similar to the measurement, 
but somewhat shifted to the upper right. This is probably due to the over-prediction of 
the maximal wave height. The RD schemes gives a better approximation of the maximal 
height, but the following graph is smoothed out stronger. Finally, in Fig. 21 we present 
solutions at the gages for different grid resolutions, computed with the FVEG scheme. 
We can clearly see that the scheme converges. The steepening of the fronts becomes more 
pronounced and the perturbations during the run-down vanish on the finer grids. How¬ 
ever, the main features are already captured on the relatively coarse grid used for the 
simulation in Fig. 20 All in all, the FVEG scheme produces a good reproduction of the 
wave, and the results are clearly competitive to other numerical schemes like the ones 
presented in [?,?]. 


6 Conclusions and Outlook 

We presented an approach to ensure positivity of the water height for general finite vol¬ 
ume schemes without affecting the global time step. This was achieved by limiting the 
outgoing fluxes of a cell whenever they would create negative water height. Physically, 
this corresponds to the absence of fluxes in the presence of vacuum. A splitting of advec- 
tive and gravity driven parts of the flux preserved the well-balancing. In the context of 
FVEG schemes, we applied these techniques to develop a positivity preserving scheme 
which is well-balanced in the presence of dry areas. The scheme can also properly handle 
sonic rarefaction waves, thanks to a new entropy correction based on the evolution oper¬ 
ators. We tested the scheme on a number of problems and in general obtained satisfying 
results. 

However, the discussion of the entropy fix (see Section [T4| revealed that in supersonic 
or transonic regimes the linearised wave cones used in the EG operator do not reflect 
the physical domain of dependence adequately. We conjecture that this is the origin of 

'*'The unstructured triangulation used in [?] consists of 19824 elements, whereas the grid used here has 18750 
cells. 
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the loss of convergence for Thacker's planar solution (see Section 5.4.2| , since here the 
velocities tangential to the boundary are larger than the (vanishing) gravitational speeds. 
Two issues should be analysed further. As mentioned above, the first is the linearisation 
strategy used in 


With the entropy fix from 3.4 we made a first step towards a 
more sophisticated strategy adapted to the state of flow. The other issue is related to 
the approximation of the resulting linearised evolution operators. The approximations 
used here and in [?] are based on the approximations from [?], where they have been 
developed for the wave equations. Now for this system the second eigenvalue is always 
zero, so the sonic cone is never shifted in space with respect to the prediction point. An 
approximation taking this shift into account should give more accurate results in the 
critical regime. 

Another possibility to improve the results is the introduction of friction terms. This 
could be helpful to control the velocities at the dry boundary by slowing down the waves 
near the shoreline. Finally we will combine the new scheme with the adaptation tech¬ 
niques presented in [?]. 
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